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ABSTRACT 


A  study  was  conducted  of  the  digital  computer  simulation 
of  the  equations  of  motion  of  the  Surface  Effect  Ship  with 
six  degrees  of  freedom  as  developed  by  the  Oceanics  Corpora- 
tion.  The  version  of  the  program  under  consideration  was  a 
simulation  of  the  Bell  Aerospace  Systems  100- B  Captured  Air 
Bubble  Craft  (CAB)  as  adapted  to  the  IBM-36O/67  digital 
computer  and  operating  in  irregular  seas.   The  objective  of 
the  study  was  to  optimize  the  simulation  by  reducing  compu- 
tation time  required  to  obtain  solutions  to  the  forces  and 
moments  acting  on  the  craft  while  maintaining  a  reasonable 
degree  of  accuracy  of  the  output  variables.   CPU  time 
savings  achieved  by  use  of  the  FORTRAN  H  compiler  are 
documented.   The  dependence  of  CPU  time  and  output  accuracy 
on  the  tolerance  levels  established  for  the  integration 
algorithm  is  discussed.   CPU  time  savings  versus  output 
accuracy  as  the  tolerance  levels  are  increased  is  presented. 
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I.   INTRODUCTION 

A.   BACKGROUND 

In  recent  years  the  U.  S.  Navy  has  developed  an  interest 
in  the  Surface  Effect  Ship  (SES)  concept  for  combatant  and 
support  type  craft.   As  a  result  of  this  interest,  Oceanics 
Incorporated  was  commissioned  by  the  U.  S.  Navy  to  develop 
a  digital  computer  simulation  of  the  six  degrees  of  freedom 
equations  for  SES  craft  that  would  yield  time  domain  outputs 
of  loads  and  motions.   This  initial  report  was  delivered  to 
the  Navy  Surface  Effects  Ships  Project  Office  (SESPO)  and 
dated  August,  1971  [Ref.  1]. 

The  simulation  program  was  designed  to  achieve  a  suffi- 
ciently accurate  and  stable  solution  to  the  loads  and  motions 
equations  for  the  modeling  of  the  craft.   The  desire  for  real- 
time solutions  was  acknowledged,  as  they  would  be  valuable 
for  manned  control  simulation  studies  and  greatly  reduce 
costly  computer  time  during  case  studies  [Ref.  1], 

The  program  was  developed  using  FORTRAN  IV  computer 
language  and  was  prepared  in  modular  form  for  ease  in  adapta- 
tion to  various  SES  configurations.   The  program  delivered 
to  the  Naval  Postgraduate  School  in  October,  1972,  was  a 
simulation  of  the  Bell  Aerospace  Systems  100-B  Captured  Air 
Bubble  Craft  (CAB),  a  100  ton  craft. 

Results  of  simulated  runs  conducted  by  Oceanics, 
utilizing  a  Control  Data  Corporation  6600  digital  computer 
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(using  the  scope  3-3  operating  system)  indicated  that 
problem  time  to  computer  time  ratios,  as  low  as  3:1  for 
calm  water  runs  and  as  large  as  1:4  for  large  amplitude, 
oblique,  irregular  sea  cases  at  high  speed,  could  be 
expected  [Ref.  2], 

The  initial  studies  at  the  Naval  Postgraduate  School, 
using  the  IBM-36O/67  computer  with  the  standard  Fortran  G 
compiler  to  solve  the  100-B  loads  and  motion  program  for 
sea  state  conditions,  resulted  in  problem  time  to  computer 
time  ratios  of  1:40  to  1:70,  depending  on  the  type  run 
simulated.   As  the  sea  state  condition  increased,  the 
computer  time  increased  to  the  point  where  for  certain  type 
runs-,  2  to  3  hours  of  computer  time  were  needed  for  problem 
solution  [Ref,  3].   For  example,  broaching  studies  simulating 
65  knots  with  sea  state  3  (irregular  waves)  at  150°  encounter 
angle  and  a  30  second  run  time,  required  121  minutes  of  CPU 
time.   It  should  be  noted  however,  that  in  the  broaching 
study,  the  output  requirements  included  a  listing  of  sidewall 
gap  each  time  subroutine  RHS  is  entered  thus  adding  a  signif- 
icant amount  of  CPU  time. 

B.   OBJECTIVES 

The  high  ratio  of  problem  time  to  computer  time  for  sea 
state  operations  demanded  that  an  in  depth  study  be  made  of 
the  Oceanics  SES  simulation  program  as  applied  to  the  IBM- 
36O/67  located  at  the  Naval  Postgraduate  School's  W.R.  Church 
Computer  Facility  with  the  following  objectives  in  mind: 


13 


1.)   Careful  analysis  of  each  subroutine  to  determine 
those  areas  where  the  bulk  of  the  computer  time  is  being 
spent. 

2.)   Determine  methods  whereby  the  computational  effi- 
ciency of  the  program,  can  be  maximized  with  minimum  loss 
in  accuracy. 

This  thesis  examines  the  SES  100-B  digital  computer 
simulation  program  on  the  IBM-36O/67  computer  with  the 
above  mentioned  objectives  in  mind  and  makes  recommendations 
for  substantially  reducing  computer  time  while  maintaining 
a  reasonable  degree  of  accuracy  with  regard  to  the  loads 
and  motion  simulation. 
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II.   INTEGRATION  METHOD 

A.   INTRODUCTION 

The  SES  simulation  program  is  based  on  a  mathematical 
description  of  the  craft  in  six  degrees  of  freedom  as 
reported  by  Kaplan,  Bentson,  and  Sargent  [Ref.  1].   The 
resulting  equations,  in  a  form  suitable  for  numerical 
integration  on  a  digital  computer,  are  given  by 


(1) 
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6  =  Q  (8) 


z  =  w 


X  =  u 


(9) 
(10) 


15 


y  =  v  (ID 

i  =  R  _  (12) 

"b  =  pa(Qin  -  Qout}  <13) 

where  X,  Y,  and  Z,  are  the  3-dimensional  spatial  coordinates; 
P,  Q,  and  R,  are  the  angular  velocities  about  the  X,  Y,  and 
Z,  axis  respectively;  I  is  the  mass  moment  of  inertia,  and 
P  is  the  force.   M,  is  the  mass  of  bubble  air  in  the  plenum; 
Qin  is  the  volumetric  quantity  of  air  flow  rate  into  the 
bubble;  Q  .  is  the  flow  rate  due  to  leakage,  and  p   is  the 
standard  atmospheric  reference  density. 

The  program  is  of  modular  design  utilizing  a  parti- 
tioning technique  on  the  craft  to  calculate  the  various 
forces  and  moments  exerted  on  the  hull  as  it  moves  through 
the  water.   The  forces  and  moments  acting  on  the  craft  are 
computed  in  various  subroutines  using  initial  conditions 
supplied  by  subroutine  INCON.   Subroutine  RHS  consists  of 
a  simple  summing  algorithm  to  compute  the_  total  of  the  forces 
and  moments  in  six  degrees  of  freedom. 

Subroutine  INTGRL  calls  subroutine  RHS  for  the  calcula- 
tion of  the  right  hand  side  of  equations  defined  by  the 
Runge-Kutta-Merson  (RKM)  algorithm  for  each  required  step 
size.   The  elements  of  each  k  vector  in  the  RKM  algorithm 
are  the  ten  variables  defined  by  equations  (1)  through  (9) 
and  equation  (13) . 
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Since  the  variables  X,  Y,  and  \\>   defined  in  equations 
(10),  (11),  and  (12),  are  expected  to  have  slow  variation, 
it  was  decided  that  a  simple  means  of  integration  using 
Simpson's  Rule  could  be  applied  to  these  equations.   These 
variables  are  calculated  in  the  main  program  [Refs.  1  and  2], 

A  complete  description  of  the  computer  program  including 
its  organization,  flow  charts,  input  and  output  description, 
program  listing,  and  a  user's  manual  for  the  operation  of 
the  program  are  available  in  Refs.  2  and  4. 

B.   SUBROUTINE  INTGRL 

The  numerical  integration  technique  used  in  subroutine 
INTGRL  is  the  Runge-Kutta-Merson  (RKM)  method  for  solution 
of  a  system  of  first  order  differential  equations  of  the 
form 

y  =  f(t,y)  (11) 

the   recursion   formula  is 


k,  +  4ka  +kR  q 

Y  ^    =  Y     +  — Jl 2.  +   0(h5)  (15) 

n+1  n  2  ■ 


where 


k,   -  |  f(tn.Yn)  (16) 


k2   =!«'„*!   .   Yn+kl>  (17) 
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k3  =|f(tn+|  ■  Yn  +  TT  +  -T   >  <l8> 


h  kn   9k. 


kH  =  If^n  +  h  >  Yn  +  "F+  V  )  (19) 

3k    9k 
k5  -  |  f(tn+h  ,  Yn  +  -ji  -  -^1  +  6k„  )        (20) 


and  h  is  the  stepsize.   The  truncation  error  in  (15)  is 
given  by 

9k.        k. 

e  =  k1  -  -^  +  Hkj,  -  -f  (21) 

Initial  conditions  are  read  in  through  subroutine  INCON 
to  establish  the  starting  point  for  each  of  the  ten  variables 
plus  a  trial  stepsize  (h)  and  a  maximum  error  tolerance  level 
for  each  integrator.   Computations  are  then  carried  out  and 
the  truncation  error  given  by  equation  (21)  is  computed.   If 
any  element  in  the  error  vector  e  is  greater  than  the  corre- 
sponding maximum  tolerance  level  then  the  trial  time  stepwise 
is  halved  and  the  computations  are  repeated.   If  the  error  is 
less  than  the  tolerance  level,  but  greater  than  1/16  this 
level,  the  computations  proceed  in  accordance  with  the  RKM 
algorithm.   If  all  elements  in  e  are  less  than  or  equal  to 
1/16  the  maximum  tolerance  level,  the  stepsize  is  doubled 
in  the  next  integration  step.   This  procedure  automatically 
adjusts  the  time  step  in  the  computation  to  reflect  the 
nature  of  the  perturbations  being  experienced  by  the  craft. 
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As  a  result,  for  slowly-varying  phenomena  (e.g.,  calm 
water,  slow  speed  simulation),  the  computations  are  carried 
out  quickly  since  larger  time  steps  are  taken,  while  high 
frequency  effects  (e.g.,  heavy  seas,  high  speed  simulations), 
will  cause  an  increase  in  computation  time. 
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III.       PROCEDURE 

A.  INTRODUCTION 

The  obvious  starting  point  for  an  efficiency  study  of 
the  SES  simulation  program  was  to  gain  a  thorough  knowledge 
of  the  force  and  moment  algorithm  used.   The  object  was  to 
study  the  interconnections  and  interactions  between  the 
various  subroutines  with  a  goal  of  arriving  at  a  more  effi- 
cient (i.e.  less  time  consuming)  solution  by  reprogramming 
and/or  better  utilization  of  the  IBM-36O  software. 

A  major  reprogramming  effort  was  considered  to  be  beyond 
the  scope  of  this  thesis.   However,  investigation  of  subrou- 
tines which  were  known  to  consume  the  bulk  of  the  computation 
time  was  conducted.   In  conjunction  with  this  investigation 
various  timer  functions  built  into  the  IBM-36O  were  utilized 
to  pin  point  time  consuming  procedures. 

B.  IBM-360  COMPILER 

The  SES  simulation  program  was  originally  adapted  to  the 
Naval  Postgraduate  School  IBM-360  computer  utilizing  the 
standard  FORTRAN  G  compiler.   This  compiler  features  a  DEBUG 
facility  and  core  optimizations  which  lends  itself  to  the 
smaller  job,  teaching  environment  prevalent  at  the  school 
[Ref.  5].   Among  other  compilers,  the  IBM-36O  computer  system 
has  available  FORTRAN  H,  an  IBM  processor  which  produces  a 
highly  optimized  code  which  is  highly  desirable  for  large 
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production  programs  such  as  the  SES  simulation  program 
[Ref.  6]. 

Compilation  of  the  program  in  FORTRAN  H  resulted  in  an 
average  CPU  time  savings  of  thirty-seven  percent  over  the 
FORTRAN  G  compiler  while  yielding  exactly  the  same  output 
variable  values.   For  example,  a  simulation  of  a  thirty 
second  run  in  a  state  three,  head  sea  compiled  in  FORTRAN  G 
required  80.5  minutes  of  CPU  time.   The  same  simulated  run 
time  compiled  in  FORTRAN  H  required  only  50.7  minutes  of 
CPU  time.   A  number  of  runs,  under  varying  initial  conditions, 
were  made  to  verify  the  CPU  time  savings  using  the  FORTRAN' 
H  compiler.   Unless  otherwise  indicated  all  further  runs 
discussed  in  this  thesis  were  made  using  the  FORTRAN  H 
compiler. 

C.   TIMER  SUBROUTINES 

The  IBM-36O/67  at  the  Naval  Postgraduate  School's  W.R. 
Church  Computer  Facility  has  several  built  in  timer  subrou- 
tines which  were  helpful  in  determining  those  areas  in  the 
program  in  which  the  bulk  of  the  CPU  time  was  being  used. 
Tv/o  such  routines  were  used  in  this  study. 

1.   Subroutine  IONUM 

Subroutine  IONUM  is  a  program  which  monitors  the 
number  of  times  an  input  or  output  device  is  accessed  during 
the  run  of  a  computer  program.   In  adapting  the  Oceanics 
simulation  program  to  the  IBM-36O,  extensive  use  of  discs 
as  temporary  storage  devices  was  made.   By  monitoring  the 
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number  of  input/output  operations  with  IONUM  a  more  effi- 
cient data  blocksize  for  the  discs  could  be  determined, 
thus  reducing  the  number  of  accesses  required. 
2.   Subroutine  PROGLOOK 

Several  areas  of  the  SES  simulation  program  were 
suspected  to  be  inefficient  and/or  requiring  excessive  time 
during  the  program  solution.   In  order  to  pin  point  these 
areas  of  the  program,  subroutine  PROGLOOK  was  used. 

PROGLOOK,  when  linked  to  the  program  under  consid- 
eration, monitors  the  sequential  flow  of  the  program  instruc- 
tions.  The  routine  "strobes"  the  instruction  flow  at  a  rate 
of  50  samples  per  second  or  approximately  once  every  4000 
instructions.   In  this  manner  a  time  history  of  the  program 
functions  was  obtained.   The  output  of  subroutine  PROGLOOK, 
which  was  of  interest  to  this  study,  is  in  the  form  of  a 
histogram  which  gives  a  plot  of  percentage  of  time  the  program 
spent  performing  various  functions.   Correlation  of  the  loca- 
tion of  these  functions  with  a  computer  map  listing  of  the 
SES  simulation  program  served  to  point  out  time  consuming 
areas  within  the  program. 

D.   INTEGRATOR  TOLERANCES 

A  prime  suspect  area  of  the  simulation  program  for  excess 
time  consumption  was  the  Runge-Kutta-Merson  (RKM)  algorithm 
used  to  solve  10  of  the  force  and  moment  differential  equa- 
tions.  If  during  a  given  time  increment,  any  one  of  the  10 
error  function  values  is  found  to  be  outside  the  given 
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tolerance,  the  stepsize  is  halved  and  all  variables  are 
computed  using  the  new  stepsize.   This  process  is  continued 
until  all  10  error  functions  are  within  their  respective 
tolerances  and  then  a  new  time  increment  is  started. 

The  above  stated  procedure  does  not  constitute  an  inord- 
inate amount  of  computer  time  under  calm  water  simulations 
as  the  variable  values  tend  to  vary  slowly  and  the  RKM 
algorithm  converges  rapidly  to  the  established  tolerance 
values.   However,  craft  perturbations  caused  by  such  condi- 
tions as  simulations  of  rapid  turns,  irregular  waves  to 
simulate  sea  state,  and  loss  of  a  propulsion  unit,  can  cause 
some  of  the  variables  to  vary  at  a  high  rate  of  change,  thus 
requiring  many  iterations  of  the  algorithm  for  convergence 
to  the  established  tolerance  levels  and  resulting  in  large 
CPU  times. 

Reference  3  established  the  standard  maximum  tolerance 
levels  for  each  integrator  according  to  procedures  set  forth 
in  Ref.  2.   These  tolerance  levels  were  considered  to  be 
the  basis  of  accuracy  desired  for  the  SES"  simulation  program 
used  at  the  Naval  Postgraduate  School  and  are  listed  in 
Table  I.   The  input  and  output  variables  for  each  of  the 
13  integrators  are  shown  in  Fig.  1. 
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Table  I 
Tolerance  Values  For  Integrators 


Integrator 
1 

Number 

Output  Variables 
U 

Tolerance  Values 

.0001 

2 

V 

.0002 

3 

W 

.000001 

k 

P 

.0000001 

5 

Q 

.0001 

6 

R 

.000001 

7 

* 

.000000001 

8  • 

e 

.000000001 

9 

z 

.0001 

10 

"b 

.0001 

1.   Steady  State  Conditions 

A  number  of  runs  were  conducted  Initially,  for  each 
simulation  study  made,  in  order  to  establish  steady  state 
conditions  for  a  given  set  of  initial  conditions.   Reasonably 
accurate  steady  state  conditions  are  essential  to  prevent 
initial  transients  in  the  simulation  program  which  can 
cause  the  program  to  abort,  or  worse  yet,  to  give  erroneous 
results.   The  key  variables  to  initialize  for  a  given  input 
are  thrust,  draft,  plenum  pressure,  and  to  a  lesser  degree 
pitch  angle. 


2H 


RKM      Runge-Kutta-Merson  Method 
SR        Simpson's  Rule 

Figure    1.      Inputs   and  Outputs   of  Integrators 
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Two  run  conditions  were  selected  to  fulfill  the 
objectives  of  this  thesis.   (1)   Sea  state  3,  40  knot  speed, 
head  seas,  and  (2)   Sea  state  3,  60  knot  speed  with  turn. 
Table  II  gives  the  initial  conditions  for  these  two  run 
conditions. 


Table 

II 

Initial  Conditions 

(Sea  State 

3) 

Speed 
(kts) 

Draft 
(In) 

Plenum 
Pressure 
(PSF) 

Thrust 
(one  Engine) 
(Lbs) 

Pitch 
Angle 
(Deg) 

40.0 

30.0 

78.0 

7000.0 

0.50 

60.0 

28.0 

80.0 

10800.0 

0.45 

2.   Integrators  3  and  10 

The  greatest  effect  of  the  sea  state  condition  on 
the  SES  program  computation  occurs  on  integrator  3  (heave 
accelerations)  and  integrator  10  (time  rate  of  change  of 
the  bubble  air  mass).   The  majority  of  the  forces  and 
moments  acting  on  the  craft  when  encountering  sea  states 
is  manifested  in  changes  in  Z  directed  accelerations  and 
plenum  volume  (i.e.,  bubble  mass). 

As  a  first  approach  to  reducing  the  CPU  time  for 
sea  state  simulations,  the  changing  of  the  implementation 
of  the  RKM  algorithm  was  considered.   It  was  felt  that  if 
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Integrator  3  and  10  could  be  Isolated  and  their  solutions 
made  to  converge  before  the  other  eight  Integrations  were 
attempted,  then  a  considerable  computation  time  savings  would 
be  gained  due  to  the  reduced  number  of  iterations  needed 
through  subroutines  INTGRL  and  RHS .   However,. due  to  the 
functional  dependence  of  heave  acceleration  and  bubble  mass 
on  some  of  the  variables,  it  was  found  that  Integrators  3 
and  10  could  not  be  made  to  converge  independently  of  the 
other  eight. 

Since  both  Integrator  3  and  10  were  extremely  sensi- 
tive to  craft  perturbations  the  decision  was  made  to  loosen 
the  tolerances  on  these  two  integrators  and  compare  the 
results  with  those  obtained  using  the  standard  tolerances 
(See  Table  I) .   It  was  felt  that  if  the  percentage  change 
in  variable  values  was  not  too  great  (less  than  10  percent), 
acceptable  results  could  be  obtained  with  a  considerable 
CPU  time  savings. 

The  two  integrator  tolerance  levels  were  increased, 
in  tandom,  by  a  factor  of  ten  in  steps  of  one  for  a  total 
of  ten  runs  under  a  given  set  of  initial  conditions.   Each 
computer  run  simulated  a  30  second  craft  run  and  parameter 
values  were  printed  out  every  0.05  seconds  for  a  sample 
rate  of  20  per  second  per  variable.   In  addition,  the  root 
mean  square  (RMS)  value  and  the  average  value  of  the 
parameters  was  calculated  over  the  entire  run  interval  for 
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comparison  with  values  obtained  from  runs  using  the  standard 
tolerances.   Using  these  results  an  error  analysis  could 
be  conducted. 

3.   Integrators  1,2  and  4-9 

Analysis  of  the  CPU  time  versus  output  accuracy 
as  observed  in  the  above  runs  indicated  an  optimum  tolerance 
of  0.000007  for  integrator  3  and  0.0007  for  integrator  10. 
It  was  then  decided  to  investigate  the  sensitivity  of  the 
other  eight  integrators.   The  method  used  for  this  investi- 
gation was  to  hold  the  above  tolerances  on  integrators  3 
and  10  while  increasing  one  of  the  other  integrator  toler- 
ances by  a  factor  of  ten  (the  other  seven  tolerances  remain- 
ing at  standard  values).   Each  integrator  was  tested  in 
succession  and  the  results  compared  to  standard  runs. 
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IV.   DISCUSSION  AND  EVALUATION  OF  RESULTS 

A.   FORTRAN  H  COMPILER 

The  efficiency  of  the  FORTRAN  H  compiler  can  best  be 
illustrated  by  the  following  example.   Reference  3  cited 
the  disparity  in  CPU  times  between  simulation  runs  using 
control  statement  00201  and  00202.   These  runs  were  made 
using  the  FORTRAN  G  compiler.   Control  statement  00201  inputs 
directly  to  subroutine  INCON  the  overall  weight  of  the  craft 
plus  the  location  of  the  center  of  gravity  and  the  mass 
moment  of  inertia  about  the  X,  Y,  Z,  and  XZ-axis.   Control 
statement  00202  inputs  distributed  weights  and  moments  and 
subroutine  INCON  must  then  compute  the  craft  weight,  CG 
and  moments.   It  would  appear  that  the  additional  computa- 
tions involved  in  using  control  statement  00202  would  require 
more  CPU  time  than  statement  00201. 

Results  of  Ref.  3  indicated  that  the  opposite  result  . 
was  obtained  and  that  CPU  time  decreased  by  approximately 
10  percent  when  control  statement  00202  was  used  as  input. 
The  reason  for  this  result  has  not  been  established  and 
should  be  looked  into  by  future  investigators.   Nonetheless, 
results  of  runs  compiled  in  FORTRAN  H  reduced  the  CPU  time 
difference  cited  in  Ref.  3  to  less  than  1  percent  which 
illustrates  the  optimization  capabilities  of  this  compiler. 


29 


B.   INTEGRATOR  TOLERANCES 

The  initial  runs  were  made  to  verify  the  sensitivity 
of  the  Runge-Kutta-Merson  algorithm  to  integrator  tolerances. 
The  runs  simulated  40  knots  speed,  state  3  seas,  head  waves, 
30  second  run,  and  confirmed  that  integrators'  3  and  10 
were  the  most  sensitive  to  craft  perturbations.   The  results 
of  these  runs  were  as  follows. 

Table  III 

CPU  Time  Comparison 

(Sea  State  3) 


Integrator 

Tolerance 

CPU  Time 

No.  3 

No.  10 

(Minutes) 

0.000001* 

0.0001* 

53.55 

0.000002 

0.0002 

38.74 

0.000003 

0.0003 

33.39 

0.000004 

0.0004 

31.38 

0.000005 

0.0005 

29.20 

0.000006 

0.0006 

28.90 

0.000007 

0.0007 

26.99 

0.000008 

0.0008 

26.19 

0.000009 

0.0009 

25.03 

0.000010 

0.0010 

24.50 

PERCENTAGE  Decrease 
In  CPU  Time 

Standard 
27.66 

37.65 
41.40 

45.47 
45.73 
49.56 
51.09 
53.26 
54.25 


*  Standard  for  comparison 
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1.   Integrator  3  and  10 

The  tolerances  on  Integrators  3  and  10  were  loosened 
during  the  above  stated  runs  while  the  other  eight  integra- 
tor tolerances  were  maintained  at  their  standard  values. 
Table  III  shows  the  results  of  these  runs  with  regard  to 
CPU  time  and  Figure  2  is  a  graphical  representation  of  the 
data. 

As  the  integrator  tolerances  were  varied  the  RMS 
and  average  values  of  the  output  variables  were  compared 
with  those  obtained  using  standard  integrator  tolerances. 
The  percentage  deviations  from  standard  values  were  computed 
in  order  to  have  a  measure  of  output  variable  accuracy  as 
a  function  of  CPU  time  savings.   Table  IV  tabulates  these 
results  for  maximum  tolerance  error  levels. 

As  was  expected,  the  highest  percentage  deviations 
occured  with  heave  acceleration  (integrator  3)  and  the 
variables  concerned  with  air  flow  (integrator  10).   All 
other  variables  had  insignificant  deviations  (less  than  1%) 
from  standard  values.   In  general,  the  percentage  deviation 
of  the  output  variables  increased  almost  linearly  as  the 
integrator  tolerances  were  loosened,  with  the  maximum 
deviation  occurring  for  maximum  tolerance  level. 

Increasing  the  Integrator  tolerance  levels  above  a 
factor  of  seven  did  not  greatly  decrease  the  CPU  time  required 
for  problem  solution  but  did  tend  to  decrease  the  accuracy 
of  the  affected  output  variables.   For  example,  the  initial 
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CPU  TIME  COMPARISON  WITH  VARYING 
INTEGRATOR  3  AND  10  TOLERANCES  AND 
PERCENTAGE  DEVIATION  OF  W 
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Table  IV 

Output  Variable  Maximum  Deviation 
Prom  Standard  Values 


Output  Variable 

Draft 

Pitch  Angle 

Plenum  Pressure 

CG  Heave  Acceleration 

Surge  Speed 

X-Displa cement 

Fan  Power 

Air  Flow  Into  Plenum 

Air  Flow  Out  of  Plenum 

Bubble  Drag 

Pitch  Rate 


Maximum  Deviations  from 
Standard  Value  (Percentage) 


RMS 

Average 

0.04 

0.07 

0.82 

0.00 

0.44 

0.29 

2.50 

0.00 

0.03 

0.02 

0.04 

0.04 

1.39 

3.33 

0.05 

0.95 

0.76 

1.22 

0.61 

0.27 

0.00 

0.00 

tolerance  increase  of  a  factor  of  seven  gives  a  CPU  time 
savings  of  26.56  minutes,  while  increasing  the  tolerance 
level  to  a  factor  of  ten  achieves  only  another  2.49  minutes 
savings  in  CPU  time.   At  the  same  time  the  percentage 
deviation  in  the  value  of  heave  acceleration  increases  from 
1.8$  to  2.5$. 
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2.   Integrators  1,2  and  4-9 

With  integrator  3  and  10  tolerances  at  0.000007  and 
0.0007  respectively,  each  of  the  other  integrator  tolerances 
were  increased  in  turn  by  a  factor  of  ten.   Simulated  run 
conditions  were  as  previously  described  and  CPU  times  were 
compared.   No  appreciable  CPU  time  savings  was  observed 
during  any  of  these  runs.   The  average  CPU  time  was  2  7.59 
minutes  for  these  runs.   The  percentage  deviation  in  the 
values  of  the  output  variable  was  negligibly  small,  indi- 
cating that  these  variables  were  much  less  susceptible  to 
the  moments  and  forces  caused  by  the  simulated  sea  state 
conditions. 

C.   DIFFERENT  TYPE  RUN  CONDITIONS 

Various  run  conditions  were  simulated  in  order  to  study 
the  effect  of  loosening  the  tolerances  on  integrator  3  and 
10  while  maintaining  standard  tolerances  on  the  other  eight 
integrators.   Each  computer  run  simulated  a  30  second  craft 
run. 

Table  V  shows  the  results  of  these  runs  with  respect  to 
CPU  time  savings. 

Comparison  of  the  results  shown  in  Table   III  and  Table 
V  indicate  that  the  percentage  of  CPU  time  savings  is  about 
the  same  regardless  of  the  type  run  simulated.   Therefore 
increasing  the  tolerance  levels  on  integrators  3  and  10 
by  a  given  amount  will  result  in  a  predictable  amount  of 
computer  CPU  time  savings. 
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Table  V 

CPU  Time  Comparison  For 

Different  Type  Run  Conditions 
(Sea  State  3) 


Speed 
(knots) 

Integrator  Tolerance 
No.   3          No.   10 

Wave  Encounter 
Angle 

CPU 

Time 

(min.) 

Percentage 

Decrease  in 
CPU  tire 

40 

0.000001*      0.0001* 

150° 

52.09 

Standard 

0.000007        0.0007 

24.73 

52.52 

0.000010        0.0010 

23.47 

54.94 

60 


60 


0.000001*   0.0001*   Partial  Turn  68.04   Standard 
0.000007   0.0007  34.34     49.24 

0.000010    0.0010  31.93     53.07 


0.000001*   0.0001*    180°  Turn 

0.000007    0.0007 
0.000010    0.0010 

*  Standard  for  comparison 


69.97   Standard 
36.36     48.03 
32.46     53.61 


In  comparing  results  of  deviations  of  output  variables 
under  different  run  conditions,  no  change  was  observed 
between  the  40  knot,  head  sea  simulation  and  the  40  knot 
150°  wave  encounter  simulation  outputs.   However,  under  the 
more  severe  conditions  of  60  knots  with  turns  applied  some 
of  the  output  variable  values  did  have  an  increased  deviation 
from  the  standard  values.   Table  VI  depicts  these  results 
for  the  variables  concerned. 
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Table  VI 

Output  Variable  Maximum  Deviation 

From  Standard  Values 
(Sea  State  3,  60  knots) 


Output  Variable   Turn  Condition 


Maximum  Deviation  From 


Standard 

Val 

ue  (Percentage 

RMS 

Average 

Plenum  Pressure 

Partial 

1.55 

0.29 

180° 

1.30 

0.05 

CG  Heave 

Partial 

4.76 

1.00 

Acceleration 

180° 

5.13 

0.00 

Fan  Power 

Partial 

0.90 

4.58 

180° 

1.41 

4.61 

Air  Flow  Out 

Partial 

0.17 

1.69 

of  Plenum 

D.   SUBROUTINE  PR0GL00K 

In  order  to  determine  other  areas  of  the  SES  simulation 
program  where  large  amounts  of  CPU  time  was  being  used. 
Subroutine  PR0GL00K  was  utilized.   Integrator  3  and  10 
tolerances  were  0.000007  and  0.0007  respectively.   A  40 
knot,  state  3  sea,  head  waves,  30  second  run  was  simulated. 
Figure  3  is  a  histogram  of  the  results  of  the  run. 

It  was  found  that  hl%   of  the  computing  time  was  spent 
in  subroutine  WAVES.   The  bulk  of  this  time  was  consumed 
in  the  wave  table.   A  surprisingly  large  amount  of  time 
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(22%)  was  spent  computing  trigonometric  function  values. 
The  bulk  of  the  time  spent  in  a  given  subroutine  was  gen- 
erally found  to  be  involved  in  repetitious  DO  LOOPs. 

The  area  labeled  "Other"  Included  various  FORTRAN  built 
in  functions  and  the  remaining  subroutines  and  functions 
of  the  SES  program  which  were  not  identified  on  the  graph. 
Each  of  these  routines  contributed  less  than  1%   of  computation 
time. 
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V.   CONCLUSION  AND  RECOMMENDATIONS 

A.   INTEGRATOR  TOLERANCES 

The  computer  run  time  dependence  of  the  SES  simulation 
program  on  the  tolerance  levels  chosen  for  the  Runge-Kutta- 
Merson  algorithm  is  clearly  demonstrated  by  the  results 
shown  in  this  thesis.   Further,  CPU  time  is  directly 
related  to  the  magnitude  of  the  tolerance  levels  of  inte- 
grators 3  and  10  and  is  relatively  independent  of  the 
tolerance  levels  established  for  the  other  eight  integrators. 

Output  variable  accuracy  as  would  be  expected  is  also 
related  to  the  tolerance  levels  of  the  integrators,  however 
the  percentage  deviations  from  standard  values  is  small  and 
the  values  obtained  could  be  considered  adequate  for  most 
purposes.   In  particular,  initial  studies  under  a  given  set 
of  conditions  would  realize  a  considerable  amount  of  CPU 
time  saving  by  using  increased  tolerance  levels,  while 
maintaining  sufficiently  accurate  solutions  to  determine 
the  relative  magnitude  of  the  forces  and  moments  acting  on 
the  craft.   As  more  precise  values  are  required,  the  tolerance 
levels  could  be  progressively  tightened  with  an  accompanying 
predictable  increase  in  CPU  time  requirements  (a  valuable 
planning  tool) . 

It  should  be  noted  that  the  individual  programmer  should 
determine  his  standard  of  accuracy  for  each  output  variable. 
Since  the  magnitude  of  the  outputs  vary  greatly,  a  small 
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percentage  deviation  in  one  value  may  be  acceptable  for  a 
given  case  study  while  being  entirely  too  large  for  another 
output  value  under  the  same  conditions.   For  example,  the 
air  flow  rate  into  the  plenum  has  a  magnitude  on  the  order 
of  5000  cubic  feet  per  second  with  a  1%   deviation  from 
standard  which  may  be  entirely  acceptable,  whereas  the  CG 
heave  acceleration  magnitude  of  0.32  "g"  units  with  1% 
deviation  may  not  be  acceptable. 

Due  to  time  limitations,  the  simulated  runs  conducted 
in  this  thesis  were  confined  to  relatively  slow  speeds  in 
the  40-60  knot  range.   The  deviations  observed  for  the 
critical  output  variables  of  heave  acceleration,  plenum 
pressure,  and  air  flow  in  and  out  of  the  plenum,  increased 
with  an  increase  in  speed  and  with  the  introduction  of  turns. 
Therefore  it  is  recommended  that  further  studies  be  made  at 
higher  speeds  and  more  severe  sea  state  conditions  to  deter- 
mine the  validity  of  results  obtained  with  increased 
tolerance  levels. 

B.   DIGITAL  PROGRAMMING  TECHNIQUES 

Time  studies  of  the  SES  Loads  and  Motion  program  using 
subroutine  PR0GL00K  indicated  that  an  appreciable  amount  of 
computer  time  was  spent  in  the  trigonometric  functions.   Most 
large  computers  use  nested  Chebyshev  Polynomials  or  Taylor 
Series  expansions  as  an  efficient  and  accurate  solution  to 
the  trigonometric  functions  [Ref.  6].   A  faster  method  of 
solution,  for  example,  would  be  a  four  place  table  lookup 
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for  the  functions.   A  careful  analysis  of  the  loss  of 
accuracy  incurred  by  this  method  would  have  to  be  conducted. 
Another  consideration  would  be  the  increased  core  require- 
ments due  to  the  addition  of  a  lengthy  table. 

It  was  noted  that  considerable  computation  time  is  spent 
in  the  execution  of  DO  LOOPs  throughout  the  program.   It  is 
recommended  that  an  in  depth  study  be  made  with  the  goal  of 
replacing  the  iterative  DO  LOOP  process  with  straight  forward 
calculations  where  possible.   An  excellent  example  of  this 
type  of  programming  is  presented  in  Appendix  A  of  this  thesis. 
A  straight  forward  programming  technique  resulted  in  the 
elimination  not  only  of  DO  LOOPs  but  a  subroutine  as  well. 

C.   HYBRID  COMPUTATIONS 

Plans  are  being  formulated  at  the  Naval  Postgraduate 
School  with  regards  to  the  feasibility  of  adapting  the  SES 
Loads  and  motion  simulation  programs  to  a  hybrid  computation 
technique.   The  development  of  linearized  equations  of 
motion  presented  in  Ref.  7,  and  the  analog  computer  simulation 
developed  in  Ref.  8,  could  serve  as  a  starting  point  for 
conversion  to  hybrid  computation. 

The  Naval  Postgraduate  School's  hybrid  computation 
facilities  include  the  XDS  9300  digital  computer  interfaced 
with  a  COMCOR  5000  analog  computer  and  an  Adage  Graphics 
Terminal.   Perhaps  the  most  formidable  problem  presented  by 
a  hybrid  conversion  would  be  the  core  storage  requirement 
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in  the  digital  computer.   The  present  storage  capabilities 
could  however  be  augmented  by  disc  and/or  tapes. 

The  parameter  scaling  problems  referred  to  in  Ref.  8 
could  be  readily  handled  in  a  hybrid  configuration.   Other 
advantages  to  hybrid  computation  of  the  Loads  and  motion 
program  include  near  real  time  solutions  and  the  versatility 
of  being  able  to  change  parameters  during  the  run.   The 
visual  representation  of  the  output  parameters  afforded 
by  the  Adage  Graphic  Terminal  would  be  a  powerful  engineering 
tool.   Hybrid  computation  would  be  particularly  useful  in 
vertical  plane  studies  such  as  heave  acceleration  analysis 
and  prediction,  and  for  control  systems  studies. 
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APPENDIX  A 

An  Example  of  Programming  Changes  to  Improve 

Efficiency 


Reference  9  presented  the  following  changes  to  the  SES 
Loads  and  Motions  Program  to  improve  the  efficiency  of  the 
program.   The  affected  subroutines  are  RHS  and  INCON  and  in 
addition  the  subroutine  DMINV  is  deleted. 

The  following  portion  of  subroutine  INCON  was  deleted: 

212  DO  211  I  =  1,6 
DO  211  N  =  1,6 

211   A(I,N)  =  0.0 

DO  213  N  =  1,3 

213  A(N,N)  =  AM 
A(4,4)  =  AIXX 
A(5,5)  =  AIYY 
A(6,6)  =  AIZZ 
A (4, 6)  =  -AIXZ 
A (6,  4)  =  -AIXZ 

AIMAX  =  AMAXI(AM, AIXX, AIYY, AIZZ, ABS(AIXZ)) 
DO  214  I  =  1,6 
DO  214  J  =  1,6 

214  A(I,J)  =  A(I,J)/AIMAX 
CALL  DMINV(A,G,D) 

DO  215  I  =  1,6 
DO  215  J  =  1,6 

215  A(I,J)  =  A(I,J)/AIMAX 
IF  (D.NE.0.0)   GO  TO  10 
WRITE  (6,216) 

STOP 


^3 


Inserted  in  place  of  the  above  programming  was  the 
following: 

AMASS I  =  1.0/AM 

D  =  1.0/(AIXX*AIZZ-AIXZ*AIXZ) 

DIXX  =  AIXX*D 

DIXZ  =  AIXZ*D 

DIZZ  =  AIZZ*D 

AIYYI  =  1.0/AIYY 

GO  TO  10 

The  IN CON  parameters  were  linked  to  subroutine  RHS  by 

the  following  COMMON  statement 

COMMON/ATRIX/AMASSI, AIYYI, DIXX, DIXZ, DIZZ 

In  subroutine  RHS  the  six  element  matrix  GF(I)  was 

deleted  and  the  following  identifiers  were  substituted  for 

the  summation  of  forces:   SUMX,  SUMY,  SUMZ,  SUMK,  SUMM, 

SUMN.   In  addition  the  following  deletion  was  made. 

DO  1  I  =  1,6 
VALUE (I)  =0.0 
DO  1  J  =  1,6 

VALUE(I)  =  VALUE(I)+A(I,J)*GF(J) 
1   CONTINUE 

Substituted  for  the  above  DO  LOOPs  was  the  following 

VALUE (1)  =  SUMX*AMASSI 

VALUE (2)  =  SUMY*AMASSI-R*U 

VALUE (3)  =  SUMZ* AMASS I 

VALUE  (il)  =  SUMK*DIZZ+SUMN*DIXZ 

VALUE (5)  =  SUMM*AIYYI 

VALUE(6)  =  SUMN*DIXX+SUMK*DIXZ 

The  indicated  changes  deleted  nine  DO  LOOPs  and  one 

subroutine  resulting  in  a  much  more  efficient  program  both 

from  a  time  and  a  core  requirement  standpoint. 
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APPENDIX  B 
Equivalent  Expressions  for  Yaw  Acceleration 

Reference  9  included  the  classic  Euler  equations  of 
motion  in  six  degrees  of  freedom  and  listed  the  assumptions 
under  which  certain  terms  of  these  equations  could  be 
disregarded  as  being  insignificant.   The  question  was 
raised  as  to  whether  the  equation  describing  Yaw  Accelera- 
tion  (R)  in  classic  Euler  terms  was  equivalent  to  the  one 
defined  by  subroutines  DMINV  and  RHS. 

This  appendix  shows  that,  under  certain  assumptions,  the 
two  expresssions  are  equivalent. 

The  classic  Euler  equations  as  presented  in  [Ref.  91  are 
as  follows. 


U  =  —  -  QW  +  RV  (1) 


F 

V  =  -£  -  RU  +  PW  (2) 

m 

P 

W  =  -5.  -  PV  +  QU  •  (3) 

m 

.  ■   F,  I           F   I  Q  I 

p  _   k   zz  n  xz zz 


J   I   -I2    II   -I2    II   -I2 

XX  ZZ    XZ       XX  ZZ    XZ       XX  zz    xz 


P(I  -I  +1)1  I 

r    xx  yy   zz'  xz   R,T  ,  ,   xz  >,  -, 

[ R(i  _  i    +  ;  j 

zz  "  zz 


(4) 
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P             PR(I         -    I      )         (R2-P2)I 
A  _      m      ,             zz          xx'                      '    xz  /crv 

y  -   _ —  +  ^ +  _ ^5j 

yy  yy  yy 

F              PQ(I         -    I       )          (P-QR)I 
R=_n_+JH_xx yy_  + _xz  (g) 

zz  zz  zz 


The  fundamental  equations  of  motion  for  six  degrees  of 
freedom  developed  in  Ref.  10  are  as  follows: 


Px   -   m[U  +  QW  -  RV  -  XQ(Q2+R2)  +  YQ(PQ-R)  +  ZQ(PR+Q) ]         (7) 


F„   =   m[V  +  RU-PW-  Y„(R2+P2)  +  ZP(QR-P)  +X„(QP+R)]         (8) 

y  (j  (a  u 


Fz    =   m[W  +  PV-QU  -ZQ(P2  +  Q2)  +XQ(RP-Q)  +  YQ(RQ+P)]         (9) 


Fk   =    IXXP+  (Izz-Fyy)QR  +  m[YG(W+PV-QU)  -  ZQ(V+RU-PW) 


+    XQYG(PR-Q)  -XGZQ(PQ+R)  + YQZG(R2-Q2)]  (10) 


Fm  =  IwQ+  (I   -I  jRP  +  m[Zp(U+QW-RV)  -XP(W+PV-QU) 
m    yy     xx  zz        u  u 

+  YQZG(QP-R)  -  YQXG(QR+P)  +XQZG(P2-R2)]         (11) 


Fn  =  IZZR+  (I  y-Ixx)PQ  +  m[XG(V+RU-PW)  -YQ(U+QW-RV) 

+  ZQXG(RQ-P)  -ZGYQ(RP+Q)  +YQXG(Q2-P2)]         (12) 


where  XP ,  YG,  and  Zr  are  distances  from  the  center  of  gravity 
to  the  origin  along  the  xs  y,  and  z  axis  respectively. 
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Because  of  symmetry  of  the  craft,  the  significant  product 
terms  of  inertia  are  the  XQZG  terms.  Under  these  conditions, 
equations  (7)-(12)  reduce  to  the  following  form. 


FY  =  m[U  +  QW-RV]  (13) 


P   =  m[V  +  RU  -  PW]  (14) 


F„  =  m[W  +  PV-  QU]  (15) 


P,  =  I   P+(I   -I   )QR-I   (PQ+R)  (16) 

k    xx     zz  yy      xzv 


F  =1   Q+(I   -I   )RP  +  I   (P2-R2)  (17) 

m    yy     xx  zz      xz 


F   =1   R+(I   -I   )PQ  +  I   (RQ-P)  (18) 

n    zz    v  yy   xx      xz 


where 


xz     G  G 


Solving  equations  (13)-(l8)  for  the  desired  derivatives 
yields : 


F 

U  =  —  -  QW  +  RV  (19) 

m 

F 
V  =  _£  _  ru  +  PW  (20) 

m 
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W  =  -£   -  PV  +  QU  (21) 
m 

P  I  —I        I 

P  =  _L_  _  (_z|__yZ)QR  +  _**  (pq+r)       _          (22) 

XX  XX            XX 

F  I   -I        I 

Q  =  J5L  -  (_2"S — ^)rp  _  =£E  (p2-r2)             (23) 

yy  yy       yy 

F  I   —I         I 

r  =  JL-  _  (_yy — ^21)pq  _   *£  (rq-p)                        (2ii) 

zz  zz         zz 


Equations  (19)-(21)  are  identical  to  equations  (l)-(3). 
Upon  substituting  equation  (24)  into  (22)  and  rearranging 
equation  (23)  it  can  be  shown  that  equations  (22)-(24)  are 
equivalent  to  equations  (4)-(6). 

If  it  is  assumed  that  P,  Q,  and  R  are  each  much  smaller 
than  one,  then  product  terms  involving  them  may  be  ignored 
and  equations  (l)-(6)  reduce  to  the  following  equations  as 
stated  in  Ref.  9- 


(25) 


• 

u 

= 

F 

X 

m 

• 

• 

V 

= 

Fy 

-w-   RU 

w 

= 

F 
z 

m 

• 

p 

= 

F.  I 
k  zz 

1  VI  -1 

XX  ZZ 

F  I 
n  xz 

2 

*        I   I 
XZ      XX  zz 

-I  2 

xz 

(26) 


(27) 


(28) 
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F 
Q   =   T2-  (29) 

yy 

F  .    I 

R  =    j2-  +   P  y&  (30) 

zz  zz 


The   equations   of  motion   as   decoded  from  subroutines 
DMINV  and  RHS   are   as    follows. 


• 

u  = 

P 

X 

m 

• 

V  = 

Fy 

•£- RU 

w  = 

F 
z 

m 

• 

p  = 

F,  I            F  I 
k  zz        ,   n  xz 

2 
I   I   -  I      II- 

XX  ZZ    XZ       XX  zz 

-i  2 

xz 

Q  = 

F 
m 

• 

R  = 

n  xx        .   k  xz 

2  2 

II        -I      c         II        -I    „ 
XX    ZZ  XZ  XX    zz  xz 


(3D 


(32) 


(33) 


(34) 


(35) 


(36) 


Equations  (25)-(29)  are  identical  to  equations  (3D-(35) 
Substituting  equation  (28)  into  equation  (30)  yields: 


^9 


A  _  n    /  k  zz  n  xz       >.   xz 

R  -  j —  +  ( ?  +  ?;  i — 

izz    I  1„„  -I  „    II   -  I      zz 
XX  ZZ    XZ      XX  zz    xz 


F     F,  I  F  I   2 

n   ,   k  xz  n  xz 

—  +  ^  +  _ _ 

zz     II    -I  I   I     -I   I   ^ 
XX  zz    xz      XX  zz     zz  xz 


p 

p.  I  ,        I 

K    xz  +    p    (    _JL_  +      xz 


I      I       -I      2         nIzz        II2-II2 
XX    ZZ  XZ  XX    zz  zz    xz 


F.I  II2-II2  +  II2 

k  XZ         .  .  -a     ,       XX  ZZ      ZZ  XZ      zz  XZ  v 

^  +  *     \    p 2 ' 

II   -  I         I   CI  I   -  I   ) 

XX  ZZ     XZ  ZZ  v  XX  ZZ     XZ 


P.  I  _          F  I 
r  =  K  xz _  +  n  xx 7  (37) 

II   -I     II   -I 

XX  ZZ    XZ      XX  zz    xz 


Equation  (37)  is  identical  to  equation  (36)  therefore,  based 
on  the  given  assumptions,  the  classical  Euler  equations 
given  in  Ref.  9  and  as  developed  from  Ref.  10  are  mathe- 
matically equivalent  to  those  defined  in  subroutines  DMINV 
and  RHS.   However  it  is  felt  that  the  validity  of  disregarding 
product  terms  should  be  more  thoroughly  investigated.   It 
is  recommended  that  equations  of  the  form  (l)-(6)  be 
incorporated  in  future  digital  simulation  studies  to  determine 
the  degree  of  contribution  of  the  cross  product  terms  on  the 
resultant  output  variables. 
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